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We consider evaluation of matrix elements with the coupled-cluster method. Such calculations 
formally involve infinite number of terms and we devise a method of partial summation (dressing) 
of the resulting series. Our formalism is built upon an expansion of the product C^C of cluster 
amplitudes C into a sum of n-body insertions. We consider two types of insertions: particle/hole 
line insertion and two-particle/two- hole random-phase-approximation-like insertion. We demon- 
strate how to "dress" these insertions and formulate iterative equations. We illustrate the dressing 
equations in the case when the cluster operator is truncated at single and double excitations. Using 
univalent systems as an example, we upgrade coupled-cluster diagrams for matrix elements with 
the dressed insertions and highlight a relation to pertinent fourth-order diagrams. We illustrate our 
formalism with relativistic calculations of hyperfine constant A{Gs) and 6si/2 electric-dipole 
transition amplitude for Cs atom. Finally, we augment the truncated coupled-cluster calculations 
with otherwise omitted fourth-order diagrams. The resulting analysis for Cs is complete through 
the fourth-order of many-body perturbation theory and reveals an important role of triple and 
disconnected quadruple excitations. 

PACS numbers: 31.15.Md, 31.15.Dv, 31.25.-v,02.70.Wz 



I. INTRODUCTION 

Coupled-cluster (CC) method is a powerful and ubiquitous technique for solving quantum many-body problem. 
Let us briefly recapitulate general features of the CC method, so we can motivate our further discussion. At the heart 
of the CC method lies the exponential ansatz for the exact many-body wavefunction 

I*.) = exp(r,)|0.) = (1 + + Irf + • • • )|0.) . (1) 

V— \ (k) (k) 

Here Tt = J^k "^i cluster operator involving amplitudes of /c-fold particle-hole excitations from the 

reference Slater determinant |0i). The parametrization |^ is derived from rigorous re-summation of many-body 
perturbation theory (MBPT) series. From solving the eigenvalue equation one determines the cluster amplitudes and 
the associated energies. While the ansatz contains an infinite number of terms due to expansion of the exponent, 

(k) 

the resulting equations for cluster amplitudes contain a finite number of terms. This simplifying property is 
unfortunately lost when the resulting wavefunctions are used in calculations of matrix elements: upon expansion of 
exponents the number of terms becomes infinite. Indeed, consider matrix elements of an operator Z, e.g., transition 
amplitude between two states 



M^J = , (2) 
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with normalization Ni = (^'^1^'^). It is clear that both the numerator and denominator have infinite numbers of 
terms, e.g., 



oo oo 




In this paper we address a question of partially summing the terms of the above expansion for matrix elements, so 
that the result subsumes an infinite number of terms. 

More specifically we are interested in transitions between states of univalent atoms, such as alkali-metal atoms. 
There has been a number of relativistic coupled-cluster calculations for these systems j3i, iJj, 5, 6, 7, 8, 9] . In particular, 
calculations H, U, IE IS ignore the non-linear terms (A > 1 and ^ > 1) in the expansion (|3J); we will designate this 
approximation as linearized coupled-cluster (LCC) method. At the same time, it is well established that for the 
univalent atoms an important chain of many-body diagrams for matrix elements comes from so-called random-phase 
approximation (RPA). A direct comparison of the RPA series and the truncated LCC expansion in Ref. [T^ leads to 
a conclusion that a fraction of the RPA chain is missed due to the omitted non-linear terms. One of the methods 
to correct for the missing RPA diagrams has been investigated in Ref. Q. These authors replaced the bare matrix 
elements with the dressed matrix elements as prescribed by the RPA method. Such a direct RPA dressing involved 
a partial subset of diagrams already included in the CC method, i.e., it leads to a double-counting of diagrams. 
To partially rectify this shortcoming, the authors of Ref. Q have manually removed certain leading-order diagrams, 
higher-order terms being doubly counted. Here we present an alternative infinite-summation scheme for RPA chain 
that avoids the double counting and thus a manual removal of the "extra" diagrams. 

In addition to the RPA-like dressing of the coupled-cluster diagrams for matrix elements, we consider another subset 
of diagrams that leads to a dressing of particle and hole lines in the CC diagrams. The leading order corrections due 
to the dressing scheme presented here arise in the fourth order of MBPT, and in this paper we present a detailed 
comparison with the relevant fourth-order diagrams. Finally, we illustrate our approach with relativistic computation 
of hyperfine-structure constants and dipole matrix elements for Cs atom. In addition to dressing corrections we 
incorporate certain classes of diagrams from the direct fourth-order MBPT calculation (as in Ref. 0)01), so that 
the result is complete through the fourth order. To the best of our knowledge, the reported calculations are the first 
calculations for Cs complete through the fourth order of MBPT. 

The paper is organized as follows. First, we present a more extensive discussion of the CC formalism in Section ITU 
Further we dress particle and hole lines in Section IIVI and discuss RPA-like dressing in Section The present 
paper may be considered as an all-order extension of forth-order calculation 10, llj, and in Section IVII we present 
an illustrative comparison with the IV-order diagrams. Finally, the designed summation schemes are illustrated 
numerically in Section IVIII and the conclusions are drawn in Section IVIIII Unless noted otherwise, atomic units, 
h = \e\ ~ nie = 1, are used throughout the paper. We follow the convention of Ref. for drawing Brueckner- 
Goldstone diagrams. 



II. COUPLED-CLUSTER FORMALISM FOR UNIVALENT SYSTEMS 



In this Section we specialize our discussion of the coupled-cluster method to the atomic systems with one valence 
electron outside closed shell core. We review various approximations and summarize the CC formalism for calculation 
of matrix elements. 

We are interested in solving the atomic many-body problem. The total Hamiltonian H is partitioned as 

H = Ho + G, (4) 

where Hq is the suitably chosen lowest-order Hamiltonian and the residual interaction G = H — Hq is treated as 
a perturbation. For systems with one valence electron outside closed-shell core, a convenient choice for Hq is the 
frozen-core (y^~'^) Hartree-Fock Hamiltonian 0|. In the following, we explicitly specify the state v of the valence 
electron, so that the proper reference eigenstate |0i) of Hq is \0y) = aj|0c), where pseudo-vacuum state |0c) specifies 
the occupied core. 

For open-shell systems a general CC parametrization reads 0] 

1*,) = {exp(T,)} 1$.) , (5) 

where curly brackets denote normal product of operators. For univalent system the above ansatz may be simplified 
to 



= exp (C) S,al\0,) = V ^ S,al\Q,) . (6) 
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Here C represents cluster operator involving (single, double, triple, etc.) excitations of core orbitals 

C = C(i) + ... = (7) 

Pma ^rn^o, H" / ^ Pmnab Ct^CL^Clb^a 4~ ' ' ' : 
ma mnab 

and Sv incorporates additional excitations from the valence state v 

5„ = 1 + 5(1) + 5(2) + . . . (8) 

m mna 

In these formulae and throughout the paper we employ the following labeling convention: indexes a,b,c,d denote 
single-particle states occupied in the core |0c) and indexes m,n,r, s,t stand for remaining (virtual/excited) orbitals. 
In this convention valence states v and w form a subset of the virtual orbitals. Finally, indexes i,j,k,l stand for 
any of the above classes of single-electron orbitals. In Eqs. (|7I8() the cluster amplitudes pij stand for single-particle 
excitations and pijki for two-particle excitations, with an apparent generalization to fc-fold excitation amplitudes. 

Dictated by the computational complexity, in most applications the cluster operator is truncated at single and 
double excitations (CCSD approximation): C « C^^) + C'^) and 5^ « 1 si^^ + si^\ A further linearized (LCCSD) 
approximation consists in neglecting non-linear terms in the expansion of exponent in Eq.ljSJ, i.e., 

|M/„)LCCSD ^ ^ ^(1) ^ ^(2) ^ ^(1) ^ ^(2)^ |o^^ _ (9) 

As discussed in the introduction, the cluster amplitudes can be found from solving a proper analog of the eigen- value 
equation. We assume that these equations are solved and in a typical application we are faced with the necessity 
of computing matrix elements, Eq. |(2Jl, between two many-body wavefunctions and l^*^). As demonstrated by 
Blundell et al. 0, so-called disconnected diagrams ^3 in the numerator and the denominator of Eq. Q cancel. Their 
final expression for the exact matrix element reads 

V i""/conn ^20) 

{[1 + (^rOconnl [1 + W)con„]}'^'' 

where the matrix element Zyj^, Eq.©, is split into core Z'^"'" and valence Z,™' contributions, the diagrams comprising 
^■corc i-jgjjjg independent of the valence indexes. The valence and core parts of the normalization factor are defined 
in a similar fashion. Notice that all the diagrams in Eq. (|10|l must be rigorously connected as emphasized by subscripts 
"conn" . Since the total angular momentum of the closed-shell core is zero, the core contribution Z^°^° vanishes for 
non-scalar (and pseudo-scalar) operators and in the following discussion we will mainly focus on Z^^. 

Blundell et al. have employed the LCCSD parametrization for the wavefunction to derive 21 diagrams for 
Z™J and 5 contributions to iV™'. The LCCSD contributions to Z""'"' can be found in Ref. Ti]. It is the goal of this 
paper to go beyond these linearized LCCSD contributions. The LCCSD approximation will provide us with "skeleton" 
diagrams that will be "dressed" due to nonlinear CC terms. We display representative LCCSD diagrams in Fig. ^ 
In a typical calculation, the dominant correction to the Hartree-Fock (HF) value arise from RPA-type diagram (a) 
and Brueckner-orbital (BO) diagrams (c) and (d). (We retain the original enumeration scheme of Ref. Q for the 
diagrams.) Here are the corresponding algebraic expressions for these LCCSD contributions 

Zwv — ^ ^ -^amPwrnva h.C.S. , (H) 



Zwv — ^ ^ ^wmPmv ^ h.C.S. 
m 

Zwv ~ ^ ^ ZjnnPmwPnv i 
mn 



where h.c.s. denotes hermitian conjugation of preceding term with a simultaneous swap of the valence indexes, w v. 
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HF (a)RPA (c)BO (d) BO+ 



FIG. 1: Dominant LCCSD contributions for the matrix elements. The double arrows represent valence state, crosses represent 
matrix elements Zij and heavy horizontal lines — cluster amplitudes. In particular, the RPA diagram involves valence doubles 
and BO diagram — valence singles. Here and below we do not draw the exchange variants for the diagrams. 



III. GENERATING OBJECT C^C 

At this point we have reviewed application of the coupled cluster method to computing properties of univalent 
systems. In the remainder of this paper we deal with mathematical object 



oo oo 




As prescribed by the Wick theorem (12l |. this expression may be simplified by contracting creation and annihilation 
operators between various parts of this expression. Very complex structures may arise, so as a preliminary construct, 
consider a product C^C. Using the Wick theorem, this product may be expanded into a sum of normal forms 

C^C = (Ctc)p + (CtC)^ + {C^C)^ + ■■■= (13) 

Co + X! {a\aj^ + ^ X! ^'3''^ \a\a]aiak^ + ■■■ 

ij ij 

Here notation (C^C)^. stands for /c-body term. The zero-body term co does not have any free particle or hole lines 
and would not contribute to connected diagrams of Z^v One-body term will lead to dressing of particle and hole 
lines, discussed in Section Hvl A part of two-body term will lead to RPA-like dressing of LCCSD diagrams for matrix 
elements, as shown in Section HVI 



IV. DRESSING PARTICLE AND HOLE LINES 



In this section we focus on one-body term of the product C^C, Eq. (|14l) . and derive all-order insertions for particle 
and hole lines. To this end it is useful to explicitly express the one-body term using particle and core labels 

{C'^C)^ = - Cbaagal + Cmnaln^^n + 
ab mn 

ma ma 

Topologically, the first term is an object where a free hole line enters some (possibly very complex) structure from 
above and another hole line leaves below. The second term has a similar structure but with particle lines. In Fig. |21 
we draw these object as rectangles with "stumps" indicating where the particle or hole line is to be attached. The 
remaining terms in Eq. H14|l have both particle and hole lines involved; we will disregard these terms in the following 
discussion. 

First we prove that given a certain CC diagram for matrix elements we may "dress" all particle and hole lines as 
shown in Fig.|2| We start with a "seed" ("bare") diagram coming from a certain set of contractions in Eq. (|12|l 

-J—{0,\a.^sUC^)^° Z (C)^°5,al;|0e)socd ■ (15) 
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FIG. 2: Schematic dressing of particle and hole lines in the CC diagrams for matrix elements. 



As a next step consider a subset of terms of Eq. H12|l . constrained a,s X ^ Xq + n, ^ = + n, n ^ 1,2,3, . . .. In these 
terms carry out contractions within a product oi Xq + n operators and fiQ + n C-operators 

Within this group there are C"^^„ x C^^+n ways to pick out n pairs of operators from the two sets, being binomial 
coefficient. Once the two strings of n operators are chosen, there are n\ possible ways to contract into pairs (C^C)^^ 
object. Finally, we contract the n resulting objects into a chain (see Fig. there are n\ possible combinations. 
Combining all these factors, we recover the original factor l/(Ao!/io!) in front of the seed diagram (|15|l . 
We may define a dressed particle-line insertion f in Fig. |21 as 



ran ran 

((^^^)p-p(^^^)p-p)p_p + ---' 



where the subscript p-p denotes that we have to keep an insertion with a single incoming and a single outgoing particle 
line, e.g., {C''C)^_^ = J2mn '^■m.nCiln^n- Noticc the absence of numerical factors in front of terms of the series; this 
fact follows from the preceding discussion. Explicitly, 

r 

This series may be generated by iteratively solving an implicit equation 

r 

The very same argument holds for dressed insertions into the hole lines. 

The derivation presented above can be generalized to include simultaneous dressing of all particle-hole lines of a 
given diagram, including the inner lines of the original "bare" object (C^C)^^ itself. Below we illustrate our dressing 
scheme in case when the cluster operator is truncated at single and double excitations. 



A. Singles-doubles approximation 



With the truncated cluster operator the hole-line insertion reads 



Cba = PmbPma + „ / ^ PmnbcPmnac (18) 



2 

cmn 



and the particle-line insertion is 

Cmn ^ / ^ PnaPn^o, ~ 7: / ^ PnrabP^nrab j (19) 



2 

abr 



where we introduced anti-symmetric quantities Pmnab — Pmnab — Pnmab- 
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Diagrammatically, 



As discussed in the first part of this Section, we may dress all the particle/hole lines according to the all-order scheme 
in Fig. 01 




I m 




^= ■ 
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a 









FIG. 3: Dressing of particle and hole lines in the singles-doubles approximation. The upper and lower panels represent dressing 
of particle and hole lines respectively. 



Algebraically, 

a abr 

^ba = Sba — ^ ^ Pmb^rna ~ ^ ^ Pmnbc^rnnac i 
m rrinc 

where we introduced dressed core cluster amplitudes 

^na = ^ ^ ^ba ^ ' Psb^sn , 
b s 

nrab - '^^rf ^^ ^ns ^ ^ca ^ Pltcd^db ■ 
t s c d 

We solve Eq. iteratively 

Smn ^ "rnn Pma^^na / , Pmrab^^nrab ' 



(20) 



(21) 



abr 



^ba ~ ^^o. ^ ^ Pmb-^raa ^ ^ Prnnbc-^rnnac: (^^) 
m mnc 

where the dressed amphtudes R^^'? are to be computed with coefficients ^f.*.'* obtained at the previous step. 

With the calculated insertions ^ we can "upgrade" the LCCSD diagrams (compare Fig. and Fig.^J. To this end 
we introduce dressed matrix elements Zij and dressed valence cluster amplitudes Rmv and Rmnva^ similar to Eq. (j21|) 



rs 

-Zab = Y.^caZ.ddb. 



(23) 



cd 
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n 



(24) 



Notice that the incoming valence Une in the valence amplitudes Pmnav and Pmv is not dressed, since it does not 
represent a free end. With these objects we may dress the LCCSD diagrams, as shown in Fig. 0] 



W 
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HF 



(a) RPA 



(c)BO 



(d) B0+ 



FIG. 4: Dressing of particle-hole lines in the dominant LCCSD diagrams, Fig. Q The circled crosses represent line-dressed 
matrix elements, Eq. II23II . and the double-lined valence cluster amplitude is the dressed amplitude, Eq. (I24II . 



Numerically, we rigorously computed the four dressed diagrams shown in Fig.^ In the remaining LCCSD diagrams, 
listed in Ref. we have replaced the bare matrix elements with the dressed matrix elements, Eq. (|23|l . Notice that 
the dressing of the Hartree-Fock diagram subsumes LCCSD diagrams 

^wv ^ ~ ^ ] PmaPwaZmv + h.C.S. , 
ma 

^wv ~ ~ ^ ^ PmnbaPnwabZmv + h.C.S. . (25) 
mnah 

from Ref. , so that these diagrams are to be discarded in the present approach. We postpone discussion of numerical 
results until Section IVlII 



V. RPA-LIKE DRESSING 



In this Section we continue with the systematic dressing of the coupled-cluster diagrams for matrix elements based 
on the topological structure of the product C^C, Eq. ((TH) . Here we focus on the two-body term of this object. The 
insertion that generates the RPA-like chain of diagrams is due to a two-particle/two-hole part of the object (C^C), 

(26) 



mna 



where we used a symmetry property Cijki = Cjuk, and Cijki = Cijki — Cijik- An analysis identical to that presented for 
line dressing in Section llVI leads us to introducing an RPA-dressed object Tnbma, 

Xibrna ^a'lalaaam^ = - ^ SmnSab ^^alalaaam^ (27) 

mnab mnab 

where subscript RPA specifies that the objects inside the brackets are to be contracted so that the result has the 
same free particle/hole ends as the original bare object {G^C!)^p^, Eq. H26|l . The sign of the leading term (with the 
Kronecker symbols) is chosen in such a way that a product of this term with a particle-hole object like ZrcOcflJ results 
in the original object. 

A detailed consideration leads to an implicit equation for the RPA-dressed particle-hole insertion 

^nbjna ^jnn^ab ^ ^ ^ncra'^rbjnc • (28) 
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kAAA/ 



FIG. 5: Graphical equation for RPA-dressed insertion. The dashed horizontal line is the dressed insertion Tntma, while the 
bare object c„bma is represented by a wavy line. 





FIG. 6: Dressing particle-hole vertex of a diagram with the RPA insertion, Fig.|3 



This equation, presented graphically in Fig. can be solved iteratively. 

The resulting insertion Tnbma may dress any particle-hole vertex of a diagram as shown in Fig. |H| As an example, 



consider dressing of particle-hole matrix elements of operator Z: 
element z„ia may be derived simply as 



The equation for the dressed matrix 



;RPA 



mbna ■^nh 7 



z^P^ - z 



rnbna^nh 1 



nb 



nb 



where we used Eq. H28f) for the dressed object %nbna- Finally we obtain a set of two equations 



z^^^ - z 

^am. — '-'am 



-RPA 

mbna bn 1 



nb 



nb 

The resulting equations resemble the traditional RPA formulae for dressed matrix elements (see, e.g., Ref. ^3), but do 
not couple z^^^ and z^^. In addition, the role of the residual Coulomb interaction of the traditional RPA formulae 
is played by the matrix elements of (C^C)pjp^ object (that is dominated by second order in the Coulomb interaction, 
see Section IvUl . In addition, we would like to emphasize that the line-dressed matrix elements, Eq. (I23|l also include 
matrix elements between core and virtual orbitals; these are to be distinguished from the RPA-like matrix elements, 
Eq. 123. 

The above derivation is only valid when dressing a single isolated vertex. More complex situation arises whenever 
a horizontal cross-cut through a diagram produces several (not just two, as in Fig. (SJ unquenched particle and hole 
lines. Then the lower particle and hole lines of the object {C'C^-^pj^ could be attached to the unquenched lines of the 
bare diagram in any order and "cross-dressing" may occur. As an illustration, consider RPA-dressing of the valence 
double contribution to the wavefunction J2mna P^^a-otn'^n'^a- It arises, for example, when the diagram Z^v (see 
Fig- da)) is cut across horizontally. The valence double contains two equivalent vertexes al^aa and aJjOa- First we 
may attach the RPA object to the m — a vertex and at the second step to the n — a vertex. Apparently, this scenario 
is not covered by Eq. 127|l . since it implies that all RPA insertions are attached to the same vertex. Nevertheless, the 
RPA-like dressing can be carried out in a straightforward fashion. To continue with the illustration, the RPA-dressed 
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valence double amplitude may be defined as follows (compare to Eq. (|27|) ) 

''^mntaO-Lanaa = ( 1 + (C^C')rpa + (('^^'^)rPA (''^^ '-^) RPa) RPA + ' ' 
mna 



val. double 



Here subscript "val. double" indicates that we select a contribution having the same free ends as the R.H.S. of the 
equation. The numerical factors in front of the individual RPA contributions are derived similarly to the line-dressing 
factors (see Section Hvj) . Explicitly, we deal with a chain of diagrams 



'^mnva — Pmnva / ^ ^nbraPmrvbi 
br 

^ ^ {pncsa^shrcpTnrvh ^ncsa^nibrcPsrvb^ "t" ' ' ' . ('^0) 

hers 

Examining the structure of the above expression we finally arrive at 

-rjRPA _ -tJRPA /o1^ 

''^mnva ~ Prnnva / ^ ^nbsa ''^msvb ; K'^-^) 
sb 

i.e., we have demonstrated how to dress the valence doubles. This example can be easily generalized to several equiv- 
alent vertexes: the Eq. (|29|l has to be rewritten using as the seed object the underlying structure of the unquenched 
lines produced by a horizontal cross-cut through a diagram. 



A. Singles-doubles approximation 

Now we specialize our discussion of the RPA-like dressing to the cluster operator truncated at single and double 
excitations. In this CCSD approximation we obtain for the bare RPA-like insertion 

Cjibma ~ ~PmbP'n,a ~ / ^ PmrbcP^rac ■ 

(32) 

rc 

Notice that one has to be careful when dealing with the first (singles x singles) term; it is represented by a disconnected 
diagram and may produce undesirable disconnected diagrams for matrix elements, Eq. I|1(J|I. 

By substituting the CCSD insertion, Ea. (|32|l . into Eq. H29|l and Eq. (|31|l . we immediately derive expressions for the 
dressed matrix elements and the RPA-dressed valence doubles. For example, 

'^mnva ~ Pmnva + ^ ^ ( ^ ^ PnracPsrbc ) T^msvb i (^3) 
sb \ rc / 

^ma ~ ^nia + ^ ^ I ^ ' PnrbcPmrac 1 ^ri6 • ("^4) 
nb \ rc / 



Here we omitted a small contribution to Cnbmai Eq. H32I) . from the product of core singles. As discussed in Section IVTI 
this contribution would arise in the higher orders (sixth order) of MBPT. 

Practically, we notice that among the dominant LCCSD diagrams shown in Fig. ^ the particle-hole vertex occurs 
only in the RPA diagram Z^} . Focusing on this particular diagram, the "upgraded" Z^} is simply 

, = E ^-™^S™.a + h-c.s. . (35) 

/ dress ■^"^^ 
tna 

Notice that the use of the RPA-dressed matrix elements to upgrade the Z^v diagram, such as 

(Z(f2) = zll^p„ + h.c.s. (36) 

\ / mcl. dress ^"^^ 
ma 
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does not lead to the identical result, because it misses dressing of the w — a vertex. Moreover, we found numerically 
that dressing of both vertexes, as in more general Eq. (|35|l . is equally important. 

Finally, it is worth noting that with the RPA-like dressing scheme proposed here, the CCSD calculations would 
recover the entire chain of RPA diagrams. However, if the calculations of the wavcfunctions are done using the 
linearized version of the coupled-cluster equations (as in Rcfs. ^^i5|,0|), a part of the RPA diagrams would still 
be missing Q- To summarize, the inclusion of the CCSD non-linear terms in the CC equations is crucial for a fully 
consistent treatment of the RPA sequence. 



VI. COMPARISON WITH THE IV-ORDER DIAGRAMS 



Coupled cluster method can be straightforwardly connected with the direct order-by-order many-body perturbation 
theory. In particular, for univalent systems, when the calculations are carried out starting from the frozen-core 
Hartree-Fock potential, the lowest-order contributions to the double excitation cluster amplitudes are 

Qninva 

Pmnva ~ j ; 

~r Syi £a 

^ 9'mnab (07\ 

Pninah ~ ~J~ ■ ' ) 

I 

Here gijki is the matrix element of the residual (beyond the HF potential) Coulomb interaction between the electrons. 
It can be shown that while the LCCSD method recovers all third-order diagrams for matrix elements, it starts missing 
diagrams in the fourth order of MBPT. Our group has investigated these 1,648 complementary IVth order diagrams 
in Refs. piilll||. Among the diagrams complementary to LCCSD matrix element contributions, there are seven terms 
(class .^1x2 {Dni) in Ref. [l^ ) due to non-linear terms in the expansion of the CCSD wavcfunctions. Namely these 
^1x2 [Dni) diagrams provide the lowest-order approximation for the dressing scheme proposed here. 

In Fig. d we explore the topological structure of the Zix2 {Dni) diagrams. All these diagrams come from various 
ways of lowest-order dressing of Z^l diagram, Eq. I|ll|) . A comparison shows that the present dressing approach 
recovers five (three line-dressed and two RPA-dressed) out of seven IVth-order diagrams. The missing diagrams are 
shown in the bottom row of Fig.|7| and we call them "stretched" and "ladder" diagrams, the names being derived from 
the structure of the highlighted dressing insertions. While the ladder diagram comes from the untreated two-body 
contribution to C^C object, Eq. H14(l . the stretched diagram involves more complex three-body contribution to C^C. 
Dressing with these two insertions can be carried out in a way similar to the line- and RPA-like dressing schemes 
discussed in this paper and is beyond the scope of the present analysis. 

We verified that in the lowest-order MBPT approximation, Eq. ^VI\ . our dressing formulae reproduce the pertinent 
IVth-order expressions, explicitly presented in Ref. [Tol |. Furthermore, in Table^we list numerical values for individual 
contributions to magnetic-dipole hyperfine structure constant (HFS) A for ground state of ^^'^Cs and El transition 
amplitude for the principal Qpi/2 ~ 6si/2 transition in Cs. Analyzing this Table we conclude that both line- and 
RPA-like dressing are equally important. Moreover, for these two particular matrix elements there is a partial 
cancelation between the line- and RPA- dressed diagrams, so that were one of the dressings omitted, the result would 
be misleading. We also observe that the untreated "ladder" diagram contributes a negligibly small fraction of the 
total. At the same time the size of the untreated "stretched" diagram indicates that (at least for the HFS constant) 
it is as important as the RPA-like and line-dressed diagrams. 

As to the numerics, the IVth-order calculations have been carried out using relativistic B-spline basis sets as 
described in Ref. 0. We used a basis set of 25 out of 30 positive-energy {si > —rrieC^) pseudo-eigenfunctions for 
each partial wave. Partial waves Si/2 — 39/2 were included in the basis. The summation over intermediate core 
orbitals was limited to eight highest energy core orbitals for the El amplitude and included all core orbitals for the 
HFS constant. The reader is referred to paper for a description of our IVth-order code. 



We verified that numerical results for dressed all-order LCCSD diagrams (see Section are consistent with the 
values for the pertinent IVth-order diagrams. As an example consider contributions to A{Qs) HFS constant. Line- 
dressing modifies the LCCSD Z^^ diagram by —4.3 MHz in a good agreement with a value of —3.8 MHz, the sum 
of the first three corresponding IVth-order diagrams from Table U Similarly, RPA-like dressing modifies the Zwi 
diagram by 4-4.4 MHz, while the sum of IVth-order RPA diagram from TableQlis 4.8 MHz. 
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core line 




part, line 




streched 



RPA- II 




ladder ^ 



FIG. 7: Topological structure of fourth-order diagrams of class Zix2 (Dni) derived in Ref. fiol . In these diagrams horizonal lines 
denote Coulomb interaction and the topological objects coming from the lowest-order approximation to C^C are highlighted. 
The upper three diagrams arise in the leading order of the line-dressing scheme and two diagrams in the middle are related to 
the RPA-dressing scheme. 



Type 


A{6s), MHz 


{6si/2||i3||6pi/2), a.u. 


core line 


-2.0 


7.0[-3] 


particle line 


-0.90 


2.5[-3] 


valence line 


-0.92 


0.2[-3] 


RPA-I 


0.005 


-9.3[-3] 


RPA-II 


4.8 


-0.05[-3] 


stretched 


-3.0 


0.1[-3] 


ladder 


0.05 


-0.008[-3] 


Total 


-1.9 


0.41[-3] 



TABLE I: Breakdown of contributions to Zix2{Dni) class of diagrams for the hyperfine-structure (HPS) constant A for the 
6si/2 state and 6si/2 — 6pi/2 electric-dipole transition amplitude for ^''"^Cs atom. The lowest-order DHF values are 1425.29 
MHz for the HPS constant and 5.278 a.u. for the transition amplitude. Notation x[y] stands for x x lO''. 



VII. NUMERICAL RESULTS AND DISCUSSION 



To reiterate discussion so far, we have developed the formalism of line- and RPA- like dressing of the coupled-cluster 
diagrams for matrix elements. Further, we reduced our general formalism to the case when the cluster operator is 
truncated at single and double excitation amplitudes. We have also verified that in the lowest order we recover the 
relevant fourth-order diagrams both analytically and numerically. In this Section we illustrate our all-order dressing 
formalism with numerical results. 

We have carried out relativistic calculations of the hyperfine-structure (HFS) constant A for the 6si/2 state and 
651/2 ~ 6pi/2 electric-dipole transition amplitude for ^^'^Cs atom. It is worth noting that matrix elements of hyperfine 
interaction and electric dipole operator allow one to access the quality of ab initio wavefunctions both close to the 
nucleus and at intermediate values of electronic coordinate. Such a test is essential for estimates of theoretical 
uncertainties of calculations of parity-nonconserving (PNC) amplitudes. The ab initio PNC amplitudes are key for 
high-accuracy probes of new physics beyond the standard model of elementary particles with atomic parity violation. 

The results of calculations are presented in Tables UTI and IIIII In these tables we augment results of the previous all- 
order calculations 6] with two types of new contributions: (i) all-order RPA- and line- dressing, outlined in Sections IIVI 
and0 and (ii) complementary fourth-order contributions, so that the results are complete through the fourth-order 



12 



perturbation theory. Also, for the HFS constant, we incorporate the most recent values of the Breit and radiative 
corrections |17i.il8i|. 



Contribution 


Value (MHz) 


DHF 


1425.29 


LCCSDpT, Coulomb 


2283.1 


Breit [17] 


4.6 


VP+SE [18] 


-9.7 


LCCSDpT Reference 


2278.0 


Dressing 




A(line-dress) 


-11.0 


A(RPA-dress) 


4.4 


Dressing total 


-6.7 


Complementary IVth-order 


Triples" 


+17.7 


ZoxsiDnl) 


-5.5 


Zix2 {Dn\), stretched 


-3.0 


^1x2 (Dni), ladder 


+0.1 


A(ZIV) total 


+9.3 


Final ab initio 


2280.6 


Experiment 


2298.2 



"The Vlth-order contributions from triple excitations are beyond those treated in the LCCSDpT approximation. 

TABLE II: Contributions to the magnetic-dipole hyperfine-structure constant A of the ground 6si/2 state of ^^^Cs. 



Contribution 


Value (a.u.) 


DHF 


5.278 


LCCSDfO] 


4.482 


LCCSDpT Reference 


4.558 


Dressing 




A (line-dress) 


0.008 


A(RPA-dress) 


-0.007 


Dressing corr. total 


-0.001 


Complementary IVth-order 


Triples" 


-0.043 


Zoxs{D„i) 


0.014 


Zix2 (Dni), stretched 


1.0[-4] 


-^1x2 (-Dni), ladder 


-8.6[-6] 


A(ZIV) total 


-0.029 


Final ab initio 


4.528 


Experiment 




Young et al. fl9] 


4.5097(45) 


Rafac et al. [20] 


4.4890(65) 


Derevianko and Porsev f21j' 


4.5064(47) 


Amiot et al. [22] 


4.5006(13) 


Amini and Gould [23]'* 


4.510(4) 



"The Vlth-order contributions from triple excitations are beyond those treated in the LCCSDpT approximation. 
'From van der Waals coefficient Cg of the ground molecular state. 
■^Photoassociation spectroscopy; this is the most accurate determinat ion. 
■^From static-dipole polarizability of 65^/2 state with method of Ref. [2^. 

TABLE III: Contributions to (651/211-0116^1/2) electric-dipole matrix element for Cs atom. 
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LCCSDpT (perturbative triples) approximation. We depart from the results of the coupled-cluster calculation 
described in Ref. These are linearized coupled-cluster calculations, with the wavefunctions truncated at single and 
double excitations from the reference Slater determinant. In addition, following Ref. Q, the perturbative effect of 
triple excitations has been incorporated into the singles-doubles equation (LCCSDpT method) . The main consequence 
of this perturbative treatment is that the resulting valence removal energies are complete through the third order 
of perturbation energies. There is also a substantial (a few per cent for Cs) improvement in the accuracy of the 
resulting LCCSDpT hyperfine constants over the LCCSD values. At the same time the theory-experiment agreement 
for the El amplitudes significantly degrades (see Table llTT|l : while the LCCSD amplitudes differ by 0.4% from 0.03%- 
accurate experimental data the more sophisticated LCCSDpT matrix elements deviate from measurements by 
as much as 1.3%. In other words, both LCCSD and LCCSDpT methods are poorly suited for calculating parity-non- 
conserving (PNC) amplitudes in -'^^'^Cs with uncertainty of a few 0.1%. It is one of the goals of this paper to establish 
a method that would provide a consistent accuracy for both HFS constants and dipole matrix elements (and thus 
PNC amphtudes). 

Procedure. First we solved the relativistic LCCSDpT equations, as described in Ref. Q. With the computed 
cluster amplitudes we calculated LCCSDpT matrix elements and recovered results published in Ref. Further, we 
solved the line-dressing equations H22|l and computed the line-dressed cluster amplitudes and matrix elements. The 
convergence rate was fast: four iterations were sufficient to stabilize the norms of the line-dressed cluster amplitudes 
at a level of a few parts per million. Finally, we solved the iterative equation for the RPA-dressed valence amplitudes, 
Ea. (|33|l . It turned out to be very computationally intensive part of the scheme and we iterated the equations only once. 

We used the computed RPA-dressed valence amplitudes in calculations of the dominant diagram Zwv , see Eq. (|35|l . 
In the remaining diagrams that involve particle-hole matrix elements, we employed RPA-dressed matrix elements 
zfj^^. A numerical iterative solution of equations for z]^^'^, Eq. H34|) required only a few iterations to converge to 
seven significant figures. 

Line dressing. In Table HVl we illustrate the importance of line dressing; in this table we present differences between 
line-dressed and bare LCCSDpT diagrams for the the hyperfine constant and El amplitude. The dressing of the 
leading order HF diagram subsumes the LCCSD diagrams Zwv and Zwl, Eq. (|25|l . A direct calculation of these 
diagrams results in Z^v + Zwl — —22.9 MHz for A{6s) and —1.75 x 10~^ a.u. for the El amplitude. These values are 
consistent with dressing-induced modifications of the HF diagram (—22.0 MHz and —1.7 x lO^'^ a.u., respectively) 

1—1 . . (a) . 

from Table IIVI The modifications of the Z^y diagram are consistent with the values of the pertinent fourth-order 
diagrams, see Section IVTI and Table |l| A large dressing correction for HFS constant comes from the diagram Z^l; it 
is nominally a fifth-order diagram. The relative importance of this diagram is not surprising since it is based upon 
Brueckner orbitals (self-energy or core polarization effect). As for the El amplitude, the line-dressing correction is 

fa) 

dominated by Zwv] i.e., it is dominated by the fourth-order contribution. For the HFS constant, a relative smallness 

fa") 

of the line-dressing correction to Zwy diagram arises due to a delicate cancelation of relatively large contributions 
from dressing of the core, particle and valence lines of the diagram (see Table QJ. Finally, in the bottom line of the 
Table HVl we present a difference between the line-dressed (all diagrams) and bare values. The dressing of the HF 
diagram plays a negligible role here, since it is dominated by the diagrams already included in the LCCSDpT values. 
The line dressing contributes at a sizable 0.5% level to the HFS constant and at 0.2% level to the El amplitude. 



Type 


yl(6s), MHz (6S1/2I 


l-Dj|6pi/2),a.u. 




-22 


-1.7[-3] 


7(a) 


0.04 


9.3[-3] 


y(c) 


-8.6 


-0.5[-3] 


7(d) 


-0.8 


-3.6[-5] 


All diaj 


p-ams — 1 1 


8.0[-3] 



TABLE IV: Line-dressing induced modifications to dominant LCCSDpT diagrams, Eq. Ulljl . for the HFS constant A for the 
6si/2 state and 6si/2 — &P\/2 El transition amplitude for ^^^Cs atom. 



RPA-like dressing. Numerically dominant contribution due to the RPA-like dressing arises for Zwi diagram, where 
we used RPA-dressed valence amplitudes. The induced correction is as large as 0.2% for both dipole amplitude and 
HFS constant. The dressing of particle-hole matrix elements (z^ zf-^^) in diagrams beyond Zwl played a relatively 
minor role, contributing at a level of only 0.01% for both test cases. 
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Complementary fourth-order diagrams. The LCCSDpT method misses certain many-body diagrams for matrix 
elements starting from the fourth-order of MBPT. These complementary corrections in the fourth order come from 
triple and disconnected quadruple (or nonlinear double) excitations. In Ref. these corrections were classified by 
the role of triples and disconnected quadruples in the matrix elements (i) an indirect effect of triples and disconnected 
quadruples on single and double excitations lumped into class Zq^s', (ii) direct contribution to matrix elements, 2^1x2! 
(iii) corrections to normalization, Znom- More refined classification reads 



Here we distinguished between valence (T^) and core (Tc) triples and introduced a similar notation for singles (S) 
and doubles {D). Notation like I?^[Tc] stands for effect of core triples (Tc) on valence doubles through an 
equation for valence doubles. The LCCSDpT method combines several diagrams from Zix2{Tv) and Zox3{Sv[Tv]) 
classes. We removed these already included diagrams from the fourth-order triples in Tables |n] and IIIII Diagrams 
Dn\ are contributions of disconnected quadruples. As discussed in Section IVII namely one of such contributions, 
■2^1x2 (Dni), provides the lowest-order approximation to our all-order dressing scheme. In Tables Hll and Hill we added 
the contributions of untreated "stretched" and "ladder" diagrams of the ^1x2 (Dni) class and also from Zq^s (Dni) 
class. The latter contribution would have been accounted for by solving the full (not linearized) CC equations. In 
our large-scale fourth-order calculations we have employed the code described in Ref. (TJ ; all the formulas for a large 
number of diagrams and the code have been generated automatically using symbolic algebra tools. While the resulting 
fourth-order corrections from triples are at the level of 1%, we notice that there are certain noticeable cancelations 
between various diagrams. Thus a complete all-order treatment of triples would be essential for attaining the next 
level of theoretical accuracy. 

Hyperfine constant results. Details of calculation of the hyperfine constant A{6si/2) are presented in Table ITU To 
clarify the role of correlations, we first incorporate a number of small but important effects into the reference value: 
Bohr-Weiskopf effect, Breit and radiative corrections. The "LCCSDpT, Coulomb" value has been computed using the 
finite nuclear size, both for determination of wavcfunctions and computing matrix elements of hyperfine interaction 
(this accounts for 0.5%.) We also include Breit corrections from Ref. [T^l; these corrections differ substantially 
from those incorporated in Ref. 4, 6| due to order-of-magnitude important correlation corrections. Finally, radiative 
corrections to magnetic-dipole hyperfine-structure constants for the ground state of alkali-metal atoms were computed 
recently by Sapirstein and Cheng |l^. They found that the vacuum polarization and self-energy (VP-I-SE) contribute 
as much as 0.4% to the ab initio value. The reader should be careful with adopting Breit values from Ref.[l5j because 
these values do not include correlation corrections (see Ref. ^3 and references therein). The final value, marked as 
"LCCSDpT Reference" deviates by 0.9% from (exact) experimental value. 

Dressing corrections partially cancel, resulting in 0.3% total dressing contribution. Fourth-order diagrams, comple- 
mentary to those already included in the LCCSDpT value are dominated by a contribution due to triple excitations 
(0.8%). We also include the "stretched" and "ladder" IV-order diagrams missed by our dressing scheme (see Sec- 
tion lVI(l . Almost all the correlation corrections are of similar sizes but of different signs, so the dressing and IV-th order 
corrections cancel, so that the final correlation correction is only 0.1%, just slightly improving the theory-experiment 
agreement when compared with the "LCCSDpT Reference" value. Our ab initio value for the HFS constant deviates 
by 0.8% from the experimental value. 

Electric- dipole 6si/2 — transition amplitude. Details of calculation for the dipole matrix element are compiled 

in Table ITTTl We do not include Breit and radiative corrections in that Table, since the Breit interaction contributes 
only 0.02% to this matrix element p^ . and radiative corrections are not known from the literature. 

There were several high-accuracy experimental determinations of (6pi/2||£'||6si/2) matrix element. We list these 
matrix elements in the bottom of Table lllll In Refs. [l^l20l| this matrix element has been extracted from the measured 
lifetime of the 6pi/2 state. Determination of Ref. is based on photoassociative spectroscopy of cold Cs atoms 
(i.e., inferred from high-accuracy measurement of molecular potentials). Another approach to extraction of dipole 
matrix elements has been proposed by us in Ref. ,21]: we exploited an enhanced sensitivity of static electric-dipole 
polarizability a{0) of the ground state and van der Waals coefficient Cg of the ground molecular state to the matrix 
elements of principal transitions (6^3/2 ||-D||6si/2} and (6pi/2||T'||6s]^/2}- Essential to the extraction of individual 
matrix elements was a high-accuracy ratio of these two dipole matrix elements measured in Ref. 12^. Based on the 
proposed method the {6pi/2\\D\\6si/2) matrix element has been deduced from high-accuracy Cq in Ref. and 
in Ref. it has been inferred from a(0) measured in that work. The most accurate matrix element comes from 
photoassociation spectroscopy |23|; their result has 0.03% accuracy and we will use that value below for calibrating 
ab initio calculations. 




(38) 
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The reference LCCSDpT El matrix element deviates from high-accuracy measurements by as much as 1.3%. The 
correlation corrections (dressing and fourth-order) computed by us improve the agreement to about 0.6% , i.e., the 
ab initio accuracy becomes comparable to that for the HFS constant. An analysis of Table ITTll shows that due to 
cancelation of line- and RPA-like-dressing corrections the overall effect of dressing is negligible for this transition 
amplitude. At the same time, the forth-order corrections due to triple excitations beyond LCCSDpT triples are very 
large, almost 1%. There corrections due to residual fourth-order RPA corrections (.^ox3(-D„i)) are also sizable, and 
tend to decrease the effect of triples. Our forth-order calculation demonstrates that a full (beyond that of LCCSDpT) 
treatment of triple excitations improves the accuracy of ab initio transition amplitudes. 

VIII. CONCLUSION 

The main two results of this work are: (i) development and application of all-order dressing formalism for matrix 
elements computed with coupled-cluster method; (ii) first calculations of matrix elements for Cs complete through 
the fourth order of many-body perturbation theory. 

To reiterate, our dressing formalism is built upon a hierarchical expansion of the product of clusters C^C into a sum 
of 71-body insertions. We considered two types of insertions: particle/hole line insertion coming from the one-body 
part of the product and two-particle/two-hole RPA-like insertion due to the two-body part. We demonstrated how 
to "dress" these insertions and formulated iterative equations. Particular attention has been paid to the singles- 
doubles truncation of the full cluster operator and we derived the dressing equations for this popular approximation. 
We have upgraded coupled-cluster diagrams for matrix elements with the dressed insertions for univalent systems 
and highlighted a relation to pertinent fourth-order diagrams. Finally, we illustrated our formalism with relativistic 
calculations for Cs atom. 

Our relativistic calculations also include a large number of fourth-order diagrams complementary to LCCSDpT 
method (Linearized Coupled-Cluster Single-Doubles method with perturbative treatment of Triples; it is the most 
sophisticated CC approximation applied in relativistic calculations for Cs so far). The resulting analysis is complete 
through the fourth-order of many-body perturbation theory. We find that these complementary diagrams substantially 
improve the theory-experiment agreement for an important electric-dipole 6si/2 — (ipi/2 transition amplitude, and 
slightly better the agreement for the hyperfine constant. We found sizable cancelations between various fourth-order 
contributions; a full all-order treatment of triple and disconnected quadruple excitations is desirable to further improve 
the theoretical accuracy. 
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